Maxwell's equations approach to soliton excitations of surface plasmonic resonances 
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We demonstrate that soliton-plasmon bound states appear naturally as propagating eigenmodes 
of nonlinear Maxwell's equations for a metal/dielectric/Kerr interface. By means of a variational 
method, we give an explicit and simplified expression for the full- vector nonlinear operator of the 
system. Soliplasmon states (propagating surface soliton-plasmon modes) can be then analytically 
calculated as eigenmodes of this non-self adjoint operator. The theoretical treatment of the system 
predicts the key features of the stationary solutions and gives physical insight to understand the 
inherent stability and dynamics observed by means of finite element numerical modeling of the 
time independent nonlinear Maxwell equations. Our results contribute with a new theory for the 
development of power-tunable photonic nanocircuits based on nonlinear plasmonic waveguides. 



Nonlinear plasmonics brings the richness of the non- 
linear dynamics into the nano scaled world of photonics. 
The interest in finding nano scaled optical solitons has 
been prominent in the last few years. Recent studies 
have demonstrated self compensated diffraction beams 
on a single metal/dielectric (MD) interface [1 j, plasmon- 
solitons in MDM configurations [2], solitons in systems 
with gain and loss [3], in waveguide arrays @H6], and in 
chains of nano particles [7] and Hydrogen atoms [8] . Typ- 
ically, the nonlinear response is directly associated to the 
strong field of the SPP close to the interface, giving rise 
to surface enhancement effects [9]. 

SPP waves represent one of the most confined form of 
light known nowadays and therefore conceived as the sig- 
nals for nano-scaled circuits [10]. Their sub- wavelength 
confinement is related to the fact that they are coupled 
light-plasma waves with a propagation constant greater 
than that of the light cone [11] and require, in general, 
complex setups for their excitation via linear [12] or non- 
linear [13] mechanisms. In the context of SPP excitation, 
we demonstrate that nonlinearities play a subtle role be- 
cause, similarly to SPP's, spatial solitons in bulk dielec- 
tric media also have a propagation constant above that of 
the plane waves at the same frequency [14, . Hence, these 
two waves can couple [15 , given that their dispersion re- 
lations intersect (see Fig. [TJd), so, even in the absence 
of strong SPP fields, a nonlinear resonant excitation of 
a SPP can occur via a monochromatic beam. In order 
to show explicitly the last statement, and without loss 
of generality, most of the SPP field is located in a linear 
dielectric layer in our analysis. 

In this letter, we analyze from first principles the prop- 
erties of the stationary states of nonlinear Maxwell's 
equations for a metal/dielectric (linear) /Kerr (MDK) in- 
terface (see Fig. flk). We use a variational ansatz, which 



allows us to diagonalize the full vector nonlinear opera- 
tor of the system and derive a simplified expression for 
it. The eigenmodes of this operator, the soliplasmons, 
can be understood as stationary nonlinear hybrid states 
formed by the coupling between a SPP and a spatial soli- 
ton (which are not orthogonal). The non-self adjoint char- 
acter, in general, of this operator, which is an inherent 
property of the vectorial nature of the system (due to the 
presence of the SPP field), has important physical impli- 
cations being responsible of an asymmetric coupling be- 
tween the plasmon and the soliton. This essential feature 
was not captured by previous heuristic models [15] intro- 
duced by analogy with coupled linear SPP's [16] further 
analyzed in the context of Josephson junctions [T7] . 

The simplified model presented here brings the relevant 
physics associated to the soliplasmon solutions, and en- 
ables the simulation of realistic devices. In particular, it 
shows (i) that the soliton tail at the metal interface exerts 
the main driving action of the coupling, (ii) anti-crossing 
in the dispersion relations, (iii) a soliplasmon resonant 
amplitude such that it becomes a nonlinear SPP wave 
as a particular solution, and (iv) a neat classification of 
the solutions in terms of the relative phase between the 
SPP and soliton components. Moreover, first principle 
modeling, based on finite element analysis, is used to 
integrate the time independent vectorial and nonlinear 
Maxwell equations in two dimensions. Simulations are 
used to analyze stability of soliplasmons and to show the 
associated dynamics, which can be understood in terms 
of the variational model. Although derived here for the 
monochromatic case, its extension to describe plasmonic 
pulses [18] will bring insight to further study soliton dy- 
namics and frequency conversion effects involving these 
hybrid plasmon-soliton waves. 

We analyze the full-vector nonlinear equation for a 
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FIG. 1. (color online), (a) metal/dielectric/Kerr (MDK) 
structure with linear dielectric constants £ m = — n^, Ed — n\ 
and sk = n 2 K , respectively, (b) Dispersion of a SPP in a 
lossy and dispersive metal (black) and a spatial soliton (gray) 
owning two different amplitudes. Circles enclose the match- 
ing points (/3 m ,cJm) and the dashed line marks the light cone 



monochromatic wave (cw) with frequency uj 
scribed by its electric field E w : 
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where C v = V [Vo] is the vector linear operator and 
Pnl = -fc 2 x (3) /3{2|E w | 2 E w +E2E*}. The MDK ge- 
ometry is illuminated from the Kerr medium with a beam 
parallel to the interface (see Fig. |TJl) and diffraction 
along y is neglected. Below, we find analytical solutions 
for this problem by means of dynamical equations derived 
from Eq.Q using a variational method. The choice of the 
variational ansatz is motivated by the heuristic model in 
Ref . [T5] . according to which soliton-SPP states are ex- 
pected to exist, and thus taken in the form: 



and the SPP behaves linearly in the absence of coupling. 
The weak coupling is automatically satisfied when the di- 
electric layer is present, although it can be also achieved 
in MK structures if a is big enough. Full-vector numer- 
ical simulations show that in this situation the field is 
mainly transverse in the nonlinear region. Thus, to a 
first approximation, u « (1,0). Our ansatz has only two 
variational parameters: c p (z) and c s (z). However, c s (z) 
appears nonlinearly in Eq. (J2|, preventing the soliplas- 
mon from behaving as a linear superposition of modes. 
Parallel illumination of the MD interface ensures paraxi- 
ality unless high energy transfer takes place. For weakly 
coupled quasi-stationary soliplasmon states this is true, 
justifying the paraxial approximation and the fact that 
the soliton position, a, is not considered as a variational 
parameter. These two approximations are supported by 
our full- vector modeling (see Figs. [4|6|. 

We substitute the ansatz Eq.Q in the paraxial version 
of Eq.Q (d 2 ->> 2ikn K d z - k 2 rr K ) for e x = xE^x, z), 
—id z e x = Me Xl and distinguish between two regions: lin- 
ear (x < d) and nonlinear (x > d), d being the dielectric 
layer width. For x < d, Pnl = and we take into ac- 
count: (i) e p is an eigenfunction of the linear operator 
with C v = — V [^^ 1 V£l°] and eigenvalue f3 2 ; (ii) (j) s is 
an eigenfunction of the scalar nonlinear operator (vector 
effects are negligible for <j> 3 ) with eigenvalue /i s , but be- 
haves linearly in this region (a > d). We then project 
the resulting equation into the plasmon component by 
multiplying it by e px and integrating over x G ]— oo,oo[. 
Since the soliton and SPP modes are not orthogonal 
(I x e pxfs 7^ 0) we obtain a coupled equation for dc p /dz 
and dc s /dz. An analogous procedure applied in the re- 
gion x > d, projecting now with respect to / s , gives a 
second equation for dc p /dz, dc s /dz. Linear combinations 
of both equations in the weak coupling approximation 
give the variational equations of the soliplasmon system 
(|c> = [c p ,c s ] T ): 



E w (z, z) = [c p (z)e p (x) + uc s (z)f s (x - a; c s (z))] e lk7lKZ . 

(2) 

The field e p (x) = [e px (x),e pz (x)] corresponds to a 
TM-SPP on the MD interface with propagation con- 
stant f3 p = k^/smSd/ [e m + £d], which is a stationary 
solution of Eq.Q with Pnl = 0, and hence the com- 
plex amplitude c p (z) = c p (0)e lflpZ (j3 p = krix + /ip). 
The dielectric function £l(x) is in this case that of 
the MD interface. The soliton term in Eq.Q, cj) s = 
Csfs = c s sech (\/knKj \c s \ [x — a]), is located at a dis- 
tance a from the MD interface and it is a solution 
of the stationary scalar (C v = 0) and paraxial non- 
linear Schrodinger equation with c s (z) = c s (0)e 2/isZ : 

jl/[2/cnx]<9^ + 7 |0s | 2 1 <t>s = /i s 0s, so its wave-number is 

Ps = kn K + ii s , with fjL 8 = 7 |c s (0)| 2 /2, 7 = kn KX i3) /2 = 
keocn2UK /2. It is considered that the overlapping of the 
SPP and soliton fields is always small (weak coupling) 
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where \i v « f3 p - kn K and fi s (\c s \) 7l c s| 2 / 2 Ws ~ 
kriK + /is = krixW + #l c s| 2 ])- Note the dynamical 
propagation constants are approximately the stationary 
ones due to the weak coupling assumption. We empha- 
size that the origin of the off diagonal terms (coupling) 
in Eq. ([3| is the non-orthogonality between the plas- 
mon and soliton. Unlike in Ref. [15 , the coupling is 
not symmetric in general (M12 7^ M21) because M is 
not a hermitian operator, in fact q = qN p /N s1 where 



q = k/[2n K N p ] j x \e L - n 2 K \ e px f s , N p = J x \e px \ 2 and 
N s ee f x fg. Eqs. Q have no free parameters, all vari- 
ational constants are given in terms of e m , Ed, £k and 
n 2 : N p = [e 2 d n d e 2 m K m ] 

where n m d 
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FIG. 2. Soliplasmon dispersion for soliton peak power levels 
(a) g\c s \ 2 = 0.2 and (b) g\c s \ 2 = 0.4. Dashed lines correspond 
to dispersion of soliton and SPP, separately (as in Fig. [T]d). 
Inset in (a) zooms in the anti-crossing region and the inter- 
section of the dotted lines marks the predicted anti-crossing 
point. Soliton distance is fixed to a = 3/ kd in both figures. 
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and is driven by the exponentially decaying soliton tail. 
As a consequence, we predict that a strong soliton in an 
MDK system will excite a weak SPP (N s > N p , \c s \ > 
\cp\) at a rate q, whereas a strong SPP will excite a weak 
soliton at the much smaller rate q ~ \c s \. 

Soliplasmons are stationary states of Eq.([3| given by 
the eigenvalues (/i) and eigenvectors of the matrix 
M. The wave number j3 = kriK + /i is determined from 



qq, 8 = 0,7r, 



(5) 



where Ji = [/i p +/i s ]/2 is the propagation constant mean 
value and A M = \/i v — /l s }/2 the soliton-plasmon wave- 
number mismatch. For an arbitrary c s o, the associated 
eigenvectors are \/jls) = c s o[q/ {/jls — /i p } , 1] T and so 
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Note from Eq. ([5| that /i$ > /i v and /i^ < so the 
plasmon term in Eq. ([6]) is positive (negative) for 8 = 
(tt). 6 has the meaning of the relative phase between the 
plasmon and soliton components. Equation ([5| provides 
the dispersion relations in the implicit form uo vs (3$ = 
ukoj/c + /is (a;; |c s q| , a), which clearly present the anti- 
crossing behavior of the J = 0, 7r branches (see Fig{2|. Far 
from the matching point {/3 m ,u; m } (see Fig. [TJd), these 
curves tend to the non-interacting soliton and plasmon 
(dotted curves) and the S = (tt) soliplasmons switch 
from soliton (plasmon) to plasmon (soliton) behavior as 
uj increases through cj m , which is tuned with the anti- 
crossing point by |c s q|- 
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FIG. 3. (color online). P vs /i/k curves for 5 = and S = tv 
soliplasmon families at uj c± 1.26 PHz (A = 1.5 /xm) for: a = 
3/im (blue), a = 4 /mi (red), and a = 5 /im (black). Insets: 
x (solid) and z (dashed) dimensionless components of typical 
mode profiles S = \J x^E for S — (right) and 5 — tt (left). 



The soliplasmon guided power is P = f x E*H y /2 = 
uso/[2f3] f x e\E x \ 2 , where e = sl + £nl- P vs /i curves 
for 0, 7r soliplasmons can be explicitly calculated using 
the variational solutions in Eqs.([5| and (|6|: 



ft 00 



CSpSd 

2n K 



| C S 1 



1 2 



/IS - 



2qq p 



J 

(7) 

where P s is the soliton power, q ps = J x e px f Sl and 
higher order overlapping terms are neglected. The 
P(/i; uj, a) curves (see Fig. [3| have the two branches /iq > 
/ip (right) and /i^ < /i p (left) separated by the vertical 
line /i = /i p . Far from this asymptote (\/i$ — /i p \ / k ^ 1) 
both branches coalesce into the monotonically increasing 
soliton curve P s (/i) ~ /i 1 ^ 2 . However, near the resonance 
/i & /ip both curves present divergent slopes, since the 
singular SPP term in Eq.([7|) dominates, and a gap in /i 
is created. According to Eq.([5|, this is proportional to 
^Jqq ~ exp (— ay/kriKj\c s \/2} so as a increases the gap 
narrows exponentially and vice versa. All these features 
are supported by the numerically computed stationary 
solutions of Eq. ([!]) (see Fig. [3|. These eigenstates of 
M are searched in the form E(x, z) = E(x)e l/32; using 
an iterative Fourier method that fixes a, so that differ- 
ent families P(/i; a) are found separately. Soliplasmons 
of 5 = 0, 7r types naturally appear (see insets in Fig. 

Stability of soliplasmons with frequency uj c± 1.26 PHz 
has been checked by numerical propagations over z = 
60 /im with an input noise of 20% in amplitude at z = 
without ohmic losses (e m G Re). The stationary states 
are computed in a single interface configuration between 
silica glass, e K = e d = 2.09, n 2 = 2.6 x 10~ 20 m 2 /W, 
and silver e m = —82. Coupling coefficients q, q are in- 
versely proportional to N p: N s and hence a weak plasmon 
(N s ^> N p ) is driven by the soliton and vice versa. We 
focus on the former case, N s ^> N p : mathematically, 
\c s \~ 1 d\c s \/dz <C \cp\~ 1 d\cp\ / dz and q <C q, so soliton 
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FIG. 4. (color online), (a) |£cc| of the n e ff = 1.464 soliplasmon 
propagating along 60 /mi. The inset on the left magnifies \8 X \ 
over the area where it is placed, showing the noise is ejected 
away. Stream lines show the power flow, which magnitude is 
represented by the arrows. Black arrows are magnified with 
respect to the white ones. Inset on the right shows the diffrac- 
tion observed in the linear propagation over the first 12/xm. 
(b) Soliton-plasmon phase, and (c) soliton and plasmon am- 
plitudes. Sharp initial jumps are due to the large input noise. 

dynamics is essentially quasi-stationary. The previous 
assumption together with c p ^ s = \c p ^ s \ exp(i0 p?s ) allow us 
to rewrite the Eq. (|3| for the SPP amplitude, \c p \, and 
the soliplasmon phase, <j) sp = <j) p — S , 

d<t>ap _ OA , \c 3 \ , d\cp\ ||., , v 
—j— « 2A M + cos0 sp , —j^ttq\c s \sm(f) sp . (8) 

The low nonlinearity of silica glass and the power lev- 
els used here (g \c s \ 2 ~ 10 -2 ) locate the soliplasmon dis- 
persion anti-crossing frequency substantially below our 
pump frequency. For our particular choice of parameters 
the stationary solutions satisfy 2A M > q \ c s \ / \c p \ > 0, so 
the right hand side of the phase equation Q is always 
positive for the initial conditions (f3 p in the MK system is 
estimated from the vertical asymptote in Figs. |3|. Eqs. 
(|8| tells us that the soliplasmon phase will initially in- 
crease regardless the type of soliplasmon and that \c p \ 
will increase (decrease) for initial 0-(tt—) soliplasmons. 
These two features are clear in all our propagation sim- 
ulations, which integrate Eq.Q with no approximation 
(see Figs. [4|6| and permit to evaluate c p and c s as the 
peak amplitudes of the plasmon and soliton component 
of the complete solution. Figure [4] shows the propagation 
of a 5 = 7r solution with n e ^ = 1.464. Apart from being 
diffraction free (compare with linear propagation inset), 
the input noise used for all the simulations introduces 
fluctuations which propagate away from the soliplasmon 
(see inset). Here the decrease of \c p \ is associated to the 
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FIG. 5. (color online). 60 /im propagation of the n e ff = 
1.472 0-soliplasmon. (a) Electric field norm |£r| 2 , (b) soliton- 
plasmon relative phase, (c) soliton and plasmon amplitudes. 

transfer of energy from the SPP to the soliton, as shown 
by the stream lines and power flux arrows. 

Propagation of a 5 = solution, see Fig. [5j shows 
a very different behavior, since the predicted initial in- 
crease of \c p \ implies that the SPP drains energy from 
the soliton. Remarkably, d\c p \/dz = provided that 
4> sp = 0,7r, as Eqs. ^ predict and shown in Fig. 
At these points we observe that the exchange of energy 
is reversed. Qualitatively, one could explain the dynam- 
ics in Fig. [5] as follows. At z = the soliton tail reaches 
the metal interface and pumps the weak SPP wave. As a 
result of this power transfer \c p \ and j3 p (n e ^^ p ) increase 
whilst \c s \ and f3 s (n e ^^ s ) decrease. This means that the 
soliton is accelerated and the SPP is slowed down, vary- 
ing their relative phase, <j) sp . At z ~ 18 /im, (j) sp = tt but 
a 7r-soliplasmon with more stable dynamics can not be 
formed due to the different velocities between the soliton 
and SPP (d(j) sp /dz 7^ 0). As they come back in phase the 
SPP returns some energy to the soliton and the initial 
parameters are approximately restored at z = 35 /im. In 
this particular example, the soliton is slowly attracted to- 
wards the interface, presumably due to the potential well 
formed by the SPP in the Kerr medium. Simulations 
suggest that S = tt soliplasmons are more stable than 
the 5 = ones, what can be associated to the initially 
smaller value of d(j) sp /dz for our conditions in the 5 = tt 
case. Indeed, Fig. [5] suggests that SPP excitation from 
an initial soliton is possible in the 2A M > q \c s \ / \c p \ > 
regime, because the soliton tail at the metal interface has 
zero relative phase with the soliton core, what results in a 
periodic transfer of energy between the soliton and SPP. 

For 2A M < q\c s \ / \c p \ < (u <C tj m ), 7r-soliplasmons 
are expected to undergo a more rapidly varying dynam- 
ics than 5 = ones. This regime requires higher sk in 
the MDK configuration and will be reported elsewhere. 
Stability here is defined from the dynamics associated to 
\c PiS \ and (j) sp along z, rather than to the rupture of the 
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FIG. 6. (color online). Soliplasmon propagation under the ini- 
tial conditions of Fig. [5] but accounting for the metal losses, 
e m -82 + i8.3. 

nonlinear states, which is not possible since both solitons 
and SPP's are stable solution separately. 

Effect of ohmic losses is shown explicitly in Fig. |6j 
which has identical initial conditions as in Fig. [5] The 
first stages of the dynamics are very similar in the two 
cases. However, now the soliton acts as a reservoir for 
the plasmon (note that |c p (30 jam)\ « |c p (0 iim)\). As 
soon as the phase reaches (j) sp = 7r, the energy transfer 
is frustrated and \c p \ drops dramatically due to losses. 
This reduces f3 p and <p sp decreases towards the initial 
value (j) sp . Exposure to metal reflection bends the soliton 
trajectory and it goes away from the interface, leaving 
behind a SPP that will be exponentially attenuated. 

In summary, we have analyzed the stationary and dy- 
namical properties of soliton-plasmon bound states, soli- 
plasmons, by means of a simplified variational model, 
obtained from Maxwell equations. This model predicts 
the relevant physics associated to these hybrid nonlinear 
waves and is in good agreement with first principle mod- 
eling. Our novel model opens the possibility to study 
new power-tunable photonic devices based on nonlinear 
soliplasmonic waveguides. This work was partially sup- 
ported by Contract No. TEC2010-15327. 
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